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Abstract 

A recent paper [J. Chem. Phys. 132 134705 (2010)] illustrated the potential of the van der Waals density functional (vdW-DF) 
method [Phys. Rev. Lett. 92, 246401 (2004)] for efficient first-principle accounts of structure and cohesion in molecular crystals. 
Since then, modifications of the original vdW-DF version (identified as vdW-DFl) has been proposed, and there is also a new 
version called vdW-DF2 [ArXiv 1003.5255], within the vdW-DF framework. Here we investigate the performance and nature of 
the modifications and the new version for the binding of a set of simple molecular crystals: hexamine, dodecahedrane, C60, and 
graphite. These extended systems provide benchmarks for computational methods dealing with sparse matter. We show that a 
previously documented enhancement of non-local correlations of vdW-DFl over an asymptotic atom-based account close to and a 
few A beyond binding separation persists in vdW-DF2. The calculation and analysis of the binding in molecular crystals requires 
appropriate computational tools. In this paper, we also present details on our real-space parallel implementation of the vdW-DF 
correlation and on the method used to generate asymptotic atom-based pair potentials based on vdW-DF. 
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1. Introduction 

Supramolecular interactions such as steric hindrance, van der 
Waals (vdW) forces, and electrostatics play a central role in 
today's biological, nano-technology, and condensed-matter re- 
search. Materials where these interactions are important can 
be categorized as sparse matter [ 1 ], since they have low elec- 
tronic density in regions essential for cohesion and response. 
The all-pervasive vdW interactions act across such low-density 
i regions. The need for robust computational tools that provide 
insight into supramolecular systems, have given impetus to the 
development of first-principle methods that properly accounts 
for the binding of sparse matter. 

Density functional theory (DFT) within the local density ap- 
proximation (LDA) or semi-local, generalized-gradient approx- 
imations (GGA) for the exchange-correlation potential has, de- 
spite its tremendous success for dense matter, failed to consis- 
tently account for the binding in sparse matter. The van der 
Waals density functional (vdW-DF) method — both in the orig- 
inal version (2J |3), termed vdW-DFl, and in a new version, 
termed vdW-DF2 (H — go beyond these local approximations 
by using a non-local functional to approximate the correlation. 
Being first-principles, it deals with the vdW forces by including 
non-local correlations from the electron response to the electro- 
dynamical field. In the vdW-DF framework the non-local cor- 
relation energy takes the form of a double- space integral: 

E?[n] = l - Jdr JjrVr)0(r,r>(r'). (1) 

The vdW-DF method inherits the excellent performance of 
GGAs for many dense matter systems, while extending the 
reach of DFT approximations to sparse matter. 




Figure 1: Molecules forming simple model crystals bound by van der Waals 
forces: hexamine, dodecahedrane, buckmeisterfullerene/C60, and graphene. 
(from top left to bottom right). 



The original version of vdW-DFl has performed well for a 
range of systems, such as binding in molecular dimers, ph- 
ysisorption on surfaces, and polymer crystals [1]. Two of us 
have shown that it accounts for the structure of molecular crys- 
tals 0. However, vdW-DFl [2] overestimates binding sepa- 
rations by 0.2 to 0.3 A. This has led several authors to sug- 
gest adjustments in the vdW-DF method. Using the S22 [6] 
data set of molecular dimers, Klimes et al \1\ optimized several 
GGA exchange flavors for use in place of the original choice of 
revPBE x [ 8 ] (the x subscript indicates the exchange-part of the 
exchange-correlation) to the correlation in vdW-DFl; here we 
test optPBE x . Cooper developed an exchange flavor (C09 x ) for 
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vdW-DFl, [9] which goes like the gradient-expansion approx- 
imation (GEA) in the slowly- varying, high-density, limit. Very 
recently, a new version vdW-DF2 was proposed [4]. It modifies 
the inner functional, that vdW-DF uses to determine the local 
plasmon frequency ~ go(r) 2 and thereby the non-local correla- 
tion. vdW-DF2 also uses a refitted version of the PW86 X ifTOl . 
This exchange functional does not induce unphysical binding 
effects in the exchange channel and is simultaneously less re- 
pulsive than revPBE x ifTTTl . 

Computational methods designed for sparse matter benefit 
from using molecular crystals as testing grounds, because of 
the wealth of accurate experimental data on crystal structures 
and lattice parameters. Molecular crystals are also model sys- 
tem for bulk sparse matter. We therefore find it pressing to 
test vdW-DF2 and the suggested modifications of vdW-DFl on 
these systems. Here, we extend our vdW-DFl study of hex- 
amine and dodecahedrane (3 with these new modifications, 
and include two important carbon allotropes: the C60 crystal 
and graphite. The latter can be viewed as a molecular crys- 
tal of graphene flakes. For the vdW-DF correlations (termed 
vdW-DFl c and vdW-DF2 c ), we find that that both vdW-DFl c 
and vdW-DF2 c are considerably enhanced over a correspond- 
ing asymptotic atom-based pair potential form at, and a few A 
beyond, typical binding separations. 

The vdW-DF calculations for molecular crystals require an 
implementation that handle periodic systems; in particular, for 
the C60 crystal, having four molecules (240 carbon atoms) per 
unit cell, a moderately efficient and parallel implementation is 
beneficial. In this paper, we discuss some details of our parallel 
code for evaluating non-local correlations, in addition to dis- 
cussing the method used to generate results in the asymptotic 
approximation ( ' app ' ) . 

2. Computational methods 

2.1. van der Waals density functional calculations 

In the vdW-DF framework, the exchange-correlation of con- 
sists of LDA correlation, GGA exchange, and the non-local cor- 
relations: 



(a) 



^vdW-DF = E U>A + £GGA + ^ 



(2) 



The total energy, £ vdW " DF , i s obtained as in Refs. <5lfT2lfT3l. 
The exchange-correlation energy is evaluated non-self consis- 
tently using the charge density, obtained in a self-consistent 
DFT calculation with the PBE d flavor of GGA (sc- 
GGA), utilizing the ultrasoft pseudo-potential plane- wave code 
DACAPO fi51 . 

The potential energy of the crystal is the difference between 
the total energy in a configuration and the energy of isolated 
molecules: E(a) = £ vdW " DF (a) - £ vdW " DF (a -> oo) . We use 
brute force to determine the optimal value of the unit cell di- 
mension (denoted a). The molecular structure, obtained in a 
sc-GGA calculation, is kept constant as a is varied. The exper- 
imental crystal symmetry specifies the molecular orientations 
[5] for all but the C60 crystal, where we use the experimental 
orientations lfT6l . 




(b) 




Figure 2: Schematics of the evaluation the six-dimensional integral with (a) and 
without (b) a radius cutoff. The box shows the unit cell, the arrow connect to 
coordinates, the filled (dashed) circle indicates a radius cutoff. 



2.2. Implementation of the non-local correlation integral 

A parallel implementation of the non-local correlation within 
the vdW-DF framework used in recent applications 112,5, 17 ] is 
described here. It evaluates the non-local correlation of Eq. ([I]) 
by using an input charge grid provided by a software package. 
Our approach shares this post-processing strategy, using code- 
independent post- GGA evaluation of the vdW-DF method, with 
the JuNoLo (18) code. 



Figure |2.2| illustrates the evaluation of the six-dimensional 
integral in real space using multiple radius cutoffs: R\ and R2 
O [T9l E01 - The inner domain defined by |r-r'| < R\ is sampled 
on a dense grid, while the outer, defined by R\ < |r - r'| < R2, 
is sampled on a grid with half the grid-point density. Since the 
integral is six-dimensional, the cpu cost is 64 times smaller |[T3l 
per volume in the inner domain than in the outer. For periodic 
systems or systems with dimensions much larger than these ra- 
diuses, the cpu cost scales linearly with number of grid points 
N 9 but with a large prefactor ~ R 3 . For tiny, non-periodic sys- 
tem, the entire volume of the system falls within the domain 
and the scaling remains N 2 . However, depending on the grid 
sampling density, a small inner cutoff R\ can be used, and the 
costly part of the integral usually scales with N. 

In discrete form, the six-dimensional integral of Eq. ([I]) can 
be written as a double sum: Ef = (AV) 2 2/ Zj n ( r d ( Pij n ( r j) » 
where the indices /, j go over the entire grid. To evaluate a sum 
in parallel, different parts of the sum are distributed on the M 
processors. A balanced load is achieved by splitting the sum 
according to summations in integer steps: 



E-5>5>Z 




proc. M 

A scheme using interpolation to express Eq. ([T]) in terms of 
convolutions and achieves ~ N log N scaling has recently been 
reported [21]. This scheme provides a significant speedup for 
medium-size systems. However, the real-space version has cer- 
tain merits. First, the total cpu cost remains nearly constant for 
any number of processors. Second, for large systems the linear 
scaling in system size N will be important. Third, the explicit 
control over the integration domain is useful both in analysis of 
binding and for systems which are nonperiodic in some or all 
spatial directions; for instance, the tertiary structure of biopoly- 
mers are characteristically non-periodic. The explicit control 



over length scales, can also be used together with coarser ac- 
counts of the vdW forces at large separations I22l[23ll20l . 

2.3. Generating asymptotic vdW-DF pair potentials using 
Bader analysis 

The asymptotic atom-based pair potential approximation for 
the vdW energy between two molecules takes the form 



-vdW-DF = V V 6 

^app J_j J_j - _ r , - 

i j r 



(4) 



where i and j are labels for the atoms in of the two molecules. 
To generate such a potential based on vdW-DF, we partition 
the charge density of the molecule among the atoms, using the 
atoms-In-Molecules (AIM) method of Bader (24), and evaluate 
the C l £ coefficients for the different atomic pairs [20 ] . This sec- 
tion details the asymptotic form of vdW-DF, the AIM method, 
and its implementation. 

The C(y coefficient for the asymptotic vdW interaction be- 
tween two fragments, A and B, takes the general form: 



ci 



3 r 

n Jo 



duaA(iu)aB(iu) . 



(5) 



In vdW-DF the polarizability of the fragment, 



„vdW-DF 



(to) 



Jd 3 r X \ m - DF (oj,r), (6) 



is obtained by integrating over the local susceptibility: 

n A (r) 



xf"- w (<o,r) 



[9q (r)2/87i] 2 - co 2 



(7) 



This form is valid both for vdW-DFl and for vdW-DF2, since 
they only differ in how they determine the inner functional spec- 
ifying the local plasmon frequency ~ qo(r) 2 (2 El [4). 

AIM partitions the the total charge-density of a molecule 
or a solid into atomic volumes Q based solely on the topo- 
logical properties of the charge-density. The surfaces separat- 
ing the atomic volumes are defined by V^(r) -1 = where 1 
is the unit vector normal to the surface. Atomic properties, 
like charge and magnetic moments 12511 , are obtained by in- 
tegrating over the individual volumes. We use this partition 
to obtain the polarizability of an atom A using Eq. ([6]), with 
«A(r) = L d 3 r'n(r')6(r-r'). The qo(r) 2 grid is generated prior 
to this partition to avoid unphysical gradients in the boundary 
regions. 

We use and extend a code implementation [26] of the algo- 
rithm proposed by Henkelman et al [27 ], to generate Bader vol- 
umes and determine the atoms-in-molecule polarizabilities of 
Eq. ([6]). As the pseudo-potential calculations generate a pseudo- 
electron density, the first step is to include the core-electron 
density to obtain the total density. The core-electron densities 
are generated upon producing the pseudopotentials. For each 
grid point, a path of steepest descent is constructed. The set of 
paths terminating at the same maximum of the charge density, 
are assigned to the same Bader volume. 



The paths are constructed as in Ref . ITTIl : First, the charge 
density of a given grid point is compared that of all 26 adjacent 
grid points. If it is larger than its neighbors, it is considered 
a maximum; if not, the algorithm proceeds to the adjacent grid 
point with the largest charge density. This procedure is repeated 
until a local maximum is reached. It also terminates if it reaches 
a grid point that belongs to a previously assigned Bader volume. 
The algorithm scales linearly with N, because it only make a 
single loop for each grid-point. 

3. Crystal structure 

The molecular crystals of this study have simple structures. 
Hexamine and dodecahedrane forms respectively a body and 
face-centered cubic (bcc, fee) with a single molecule per unit 
cell [5]. At low temperatures, the C60 crystal is a simple cu- 
bic, with four molecule per cell. At higher temperatures, the 
molecules rotate freely and the crystal becomes effectively an 
fee |[T6ll . Graphite is a layered material, with two graphene 
sheets stacked in an AB pattern. 

Figure[3]shows the binding curve of the C60 crystal as a func- 
tion of the lattice parameter a of the simple cubic. For this 
crystal, vdW-DFl binds stronger than vdW-DF2, but at a larger 
lattice constant. The failure of DFT within GGA illustrate the 
need for non-local correlations to account for the binding in 
these crystals. 
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Figure 3: The binding curve of the C60 crystal. The lower [upper] (dashed) 
curve gives the vdW-DFl [PBE] (vdW-DF2) result. The insert shows the simple 
cubic unit cell of the crystal, which becomes face-centered cubic only if internal 
orientation are neglected. 

Table [3] shows the lattice constants and cohesion energies ob- 
tained using the vdW-DF method. The lattice constants pre- 
dicted by vdW-DFl overestimates the experimental values by 
about 0.3 A (per sheet for graphite), which is consistent with 
earlier studies [1]. The cohesion energies compares well with 
experimental values. The use of C09 x and optPBE x as exchange 
partner to the vdW-DFl c improves lattice constants. The former 



Table 1 : vdW-DF predictions of lattice parameters and cohesion energy the crystals of hexamine, dodecahedrane, C60, and graphite, compared with experimental 
values, for different combination of exchange and correlation. The experimental lattice parameters are based on low temperature measurements, except for dodec- 
ahedrane. The cohesion energies of the cage crystals (graphite) are given per molecule (atom). The bold letters identifies the results for vdW-DFl and vdW-DF2 
versions. 



Functional / molecule: Hexamine 




Dodecahedrane 


Buckyball (C60) 


Graphite 




Correlation 


Exchange a (A) 


EcohieV) 


a (A) 


EcohieY) 


a (A) 


E C oh(eV) 


c(A) 


EcohfcV) 


vdW-DFl c 


revPBE x 7.16 


1.01 




10.92 


1.46 


14.38 


1.70 


7.24 


0.053 




C09 x 6.92 


1.42 




10.44 


1.65 


14.10 


1.98 


6.44 


0.073 




optPBE x 6.96 


1.21 




10.64 


1.75 


14.22 


2.02 


6.88 


0.064 


vdW-DF2 c 


PW86 X 6.96 


0.93 




10.64 


1.35 


14.30 


1.30 


6.96 


0.053 


Exp. 


6.91 l 


0.83 2 




10.60 3 


- 


14.04 4 


1.6-1.9 5 


6.67 6 


0.052 7 


^ef. (20, 2 Ref. (29), 3 Ref. (30), 5 Ref. (31], 4 Ref. 


|16|, 6 Ref. 


(32), 


7 Ref. |33), 













is almost spot-on, but slightly underestimates them, while the 
latter overestimates them. Since the experimental value for Do- 
decaherane is based on room-temperature results, C09 x com- 
pares well with experiments also for this molecule. In both 
these modifications of vdW-DFl, the predicted cohesion ener- 
gies are quite large compared to the experimental ones. vdW- 
DF2 improves the lattice constants over vdW-DFl, giving an 
overestimation similar to optPBE x , while the cohesion energies 
are somewhat reduced. 

vdW-DF2 has the the best overall performance. The calcu- 
lated vdW-DF2 cohesion energy for the C60 crystal is lower 
than experimental observations. For this crystal, internal ori- 
entations are not optimized and this might contribute to the in- 
creased difference with the experimental data for all the mod- 
ifications. The C60 molecule also differs from the others by 
having a low surface to volume ratio. 

4. Enhancement over asymptotic pair-potential form 

Previous studies (5j [221 E21 have shown that the non-local 
correlation of vdW-DFl is enhanced compared to an atomic- 
based asymptotic pair-potential ('app') account (Eq.[4]) at bind- 
ing separation and a few A beyond. By construction 'app' does 
not include the image-plane and multi-pole effects inherent in 
the density-functional framework of vdW-DF. The importance 
of image planes were also discussed in the seminal work of 
Zaremba and Kohn in a surface-physics context |34|. Asymp- 
totic atom-based approximations are often used in force-field 
methods and in semi-empirical methods that add these interac- 
tions on top of GGA based DFT calculations (DFT-D). Since 
vdW-DF2 modifies the account of non-local correlation, we 
here investigate how this affects the enhancement over 'app'. 

Figure [4] shows the ratio between the non-local correlation 
and its corresponding 'app', AEf/E^ pp , for both vdW-DFl and 
vdW-DF2. Because of residual noise in our determination of 
the Ce coefficients, we have (< 10 %) adjusted the curves to 
reach unity at d > 16 A. The full (dashed) line gives the result 
for vdW-DFl (vdW-DF2). Both version exhibit a significant 
enhancement, with the vdW-DF2 result being somewhat larger 
and shifted to smaller separations. The less than unity ratio at 
separations about 1A shorter than the binding separations re- 
flects the built-in damping of the vdW forces in the vdW-DF 



framework; no ad hoc damping parameter needs to be intro- 
duced in the vdW-DF method. 




vdW-DF1 
vdW-DF2 



Figure 4: Role of of image-plane and multi-pole effects in the binding of a 
dodecahedrane dimer: The enhancement of non-local correlation of vdW-DF (1 
and 2) over the corresponding asymptotic atom-based pair-potential form. The 
full (dashed) line give the result for vdW-DFl (vdW-DF2). The vertical lines 
indicate the nearest-neighbor binding separation in the dodecahedrane crystal. 



5. Discussion and summary 

The performance of variations of the vdW-DF method has 
been investigated for several molecular crystals. The modifi- 
cations of vdW-DFl, involving only the exchange functional 
(C09 x , optPBE x ), improves the lattice constants, but cohesion 
energies seem to worsen for molecular crystals. That optPBE x 
did not outperform the other vdW-DFl combinations indicate 
that fitting to S22, by itself, does not guarantee high precision 
for other types of systems. C09 x gives excellent lattice parame- 
ters. The good performance of vdW-DF2 is encouraging, point- 
ing to higher accuracy for the vdW-DF method. 

The enhancement of non-local correlation over 'app' shows 
that the nature of binding both for vdW-DFl and vdW-DF2 dif- 
fers from the asymptotic 'app' description. This result further 
questions the often used 1/r 6 atom-based approximation for van 
der Waals forces when used within a few A of the binding sep- 
arations 1511221. 
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